import pandas as pd
import numpy as np
# Visualisation libraries
## Test
from IPython.display import display, Markdown, Latex
## plotly
from plotly.offline import init_notebook_mode, iplot
import plotly.graph_objs as go
import plotly.offline as py
from plotly.subplots import make_subplots
import plotly.express as px
import warnings
warnings.filterwarnings('ignore')
In these methods, a smaller interval that contains a root is identified iteratively. As these methods do not require any derivatives in general, they are also known as no derivative methods. Some of examples for these methods are:
This method is used for estimatiion of the solution of the equation $f(x) = 0$ for the real variable $x$. Here, $f$ is required to be a continuous function defined on a closed interval $[a, b]$ and $f(a)$ and $f(b)$ have opposite signs. In this case, by the intermediate value theorem, $f$ must have at least one root in the interval $(a, b)$ (for further details see [1-2]).
def bisection(f, a, b, Nmax = 1000, TOL = 1e-4, Frame = True):
# endpoint values a, b,
# tolerance TOL,
# maximum iterations NMAX
if f(a)*f(b) >= 0:
print("Bisection method is inapplicable .")
return None
# let c_n be a point in (a_n, b_n)
Nmax=Nmax+1
an=np.zeros(Nmax, dtype=float)
bn=np.zeros(Nmax, dtype=float)
cn=np.zeros(Nmax, dtype=float)
fcn=np.zeros(Nmax, dtype=float)
# initial values
an[0]=a
bn[0]=b
for n in range(0,Nmax-1):
cn[n]=(an[n] + bn[n])/2
fcn[n]=f(cn[n])
if f(an[n])*fcn[n] < 0:
an[n+1]=an[n]
bn[n+1]=cn[n]
elif f(bn[n])*fcn[n] < 0:
an[n+1]=cn[n]
bn[n+1]=bn[n]
else:
print("Bisection method fails.")
return None
if (abs(fcn[n]) < TOL):
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
n=n+1
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
Example: Finding a root of the following polynomial. $$f(x)=x^{3}-x-2.$$
First, two numbers $a$ and $b$ have to be found such that $f(a)$ and $f(b)$ have opposite signs. Given
# defining function f using Lambdas
f = lambda x: x**3 - x - 2
For the above function,
a=1
display(Latex(r'$(a, f(a))$ = (%.2f, %.2f)' % (a,f(a))))
b=2
display(Latex(r'$(b, f(b))$ = (%.2f, %.2f)' % (b,f(b))))
Now since $f(a)<0$ and $f(b)>0$, we can implement the method. Let $N_{max}=20$ and the tolerance be $10^{-4}$
def highlight_min(s):
s = np.abs(s)
is_min = s == s.min()
return ['background-color: SpringGreen' if v else '' for v in is_min]
def Method_Plot(data, Name, YL, CL = 'MediumBlue'):
display(data.style.format({'fcn': "{:.4e}"}).apply(highlight_min, subset=['fcn']))
fig = go.Figure()
fig.add_trace(go.Scatter(x = data.index, y = data.fcn, mode='lines+markers', marker_color= CL))
fig['layout']['xaxis'].update(range=[data.index[0]-0.08, data.index[-1]+0.08])
fig['layout']['yaxis'].update(range=YL)
fig.update_xaxes(showline=True, linewidth=1, linecolor='Lightgray', mirror=True,
showgrid=True, gridwidth=1, gridcolor='Lightgray', title_text = '$n$')
fig.update_yaxes(showline=True, linewidth=1, linecolor='Lightgray', mirror=True,
showgrid=True, gridwidth=1, gridcolor='Lightgray', title_text = '$f(c_n)$')
fig.update_xaxes(zeroline=True, zerolinewidth=1, zerolinecolor='Lightgray')
fig.update_yaxes(zeroline=True, zerolinewidth=1, zerolinecolor='Lightgray')
fig.update_layout(legend_orientation='v', plot_bgcolor= 'white', height= 600, width= 700)
fig.update_traces(marker_line_color= 'Black', marker_line_width=0.5, opacity=1)
fig.update_layout(legend=dict(x=.02, y=.98, font=dict(size=12, color="black"), bordercolor="Blue", borderwidth=1))
fig.update_layout(title={'text': '<b>' + Name + '<b>', 'x':0.5, 'y':0.90, 'xanchor': 'center', 'yanchor': 'top'})
fig.show()
data = bisection(f, a, b, Nmax = 20, TOL = 1e-4)
#
Method_Plot(data, Name = 'Bisection method', YL = [-.5, 2])
| an | bn | cn | fcn | |
|---|---|---|---|---|
| 0 | 1.000000 | 2.000000 | 1.500000 | -1.2500e-01 |
| 1 | 1.500000 | 2.000000 | 1.750000 | 1.6094e+00 |
| 2 | 1.500000 | 1.750000 | 1.625000 | 6.6602e-01 |
| 3 | 1.500000 | 1.625000 | 1.562500 | 2.5220e-01 |
| 4 | 1.500000 | 1.562500 | 1.531250 | 5.9113e-02 |
| 5 | 1.500000 | 1.531250 | 1.515625 | -3.4054e-02 |
| 6 | 1.515625 | 1.531250 | 1.523438 | 1.2250e-02 |
| 7 | 1.515625 | 1.523438 | 1.519531 | -1.0971e-02 |
| 8 | 1.519531 | 1.523438 | 1.521484 | 6.2218e-04 |
| 9 | 1.519531 | 1.521484 | 1.520508 | -5.1789e-03 |
| 10 | 1.520508 | 1.521484 | 1.520996 | -2.2794e-03 |
| 11 | 1.520996 | 1.521484 | 1.521240 | -8.2891e-04 |
| 12 | 1.521240 | 1.521484 | 1.521362 | -1.0343e-04 |
| 13 | 1.521362 | 1.521484 | 1.521423 | 2.5935e-04 |
| 14 | 1.521362 | 1.521423 | 1.521393 | 7.7956e-05 |
Note that the difference between $x_n$ and a solution $x$ is bounded by
$$|c_{n}-c|\leq {\frac {|b-a|}{2^{n}}}.$$Therefore, The number $N_{max}$ of iterations needed to achieve a required tolerance $\epsilon$ is bounded by
$$N_{max}\leq \left\lceil \log _{2}\left({\frac {|b-a|}{\epsilon }}\right)\right\rceil ,$$def bisection_nmax(a, b, TOL): return int(np.ceil(np.log2(np.abs(b-a)/TOL)))
In our example, we need at least iterations.
bisection_nmax(a, b, TOL = 1e-4)
14
Thus, we can modify our algorithm and remove Nmax from it.
def bisection_mod(f, a, b, TOL = 1e-4, Frame = True):
# endpoint values a, b,
# tolerance TOL,
if f(a)*f(b) >= 0:
print("Bisection method is inapplicable .")
return None
# let c_n be a point in (a_n, b_n)
Nmax= int(np.ceil(np.log2(np.abs(b-a)/TOL)))+2
an=np.zeros(Nmax, dtype=float)
bn=np.zeros(Nmax, dtype=float)
cn=np.zeros(Nmax, dtype=float)
fcn=np.zeros(Nmax, dtype=float)
# initial values
an[0]=a
bn[0]=b
for n in range(0, Nmax-1):
cn[n]=(an[n] + bn[n])/2
fcn[n]=f(cn[n])
if f(an[n])*fcn[n] < 0:
an[n+1]=an[n]
bn[n+1]=cn[n]
elif f(bn[n])*fcn[n] < 0:
an[n+1]=cn[n]
bn[n+1]=bn[n]
else:
print("Bisection method fails.")
return None
if (abs(fcn[n]) < TOL):
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
n=n+1
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
Testing the above example,
data = bisection_mod(f, a, b, TOL = 1e-4)
#
Method_Plot(data, Name = 'Bisection method', YL = [-.5, 2], CL = 'DarkGreen')
| an | bn | cn | fcn | |
|---|---|---|---|---|
| 0 | 1.000000 | 2.000000 | 1.500000 | -1.2500e-01 |
| 1 | 1.500000 | 2.000000 | 1.750000 | 1.6094e+00 |
| 2 | 1.500000 | 1.750000 | 1.625000 | 6.6602e-01 |
| 3 | 1.500000 | 1.625000 | 1.562500 | 2.5220e-01 |
| 4 | 1.500000 | 1.562500 | 1.531250 | 5.9113e-02 |
| 5 | 1.500000 | 1.531250 | 1.515625 | -3.4054e-02 |
| 6 | 1.515625 | 1.531250 | 1.523438 | 1.2250e-02 |
| 7 | 1.515625 | 1.523438 | 1.519531 | -1.0971e-02 |
| 8 | 1.519531 | 1.523438 | 1.521484 | 6.2218e-04 |
| 9 | 1.519531 | 1.521484 | 1.520508 | -5.1789e-03 |
| 10 | 1.520508 | 1.521484 | 1.520996 | -2.2794e-03 |
| 11 | 1.520996 | 1.521484 | 1.521240 | -8.2891e-04 |
| 12 | 1.521240 | 1.521484 | 1.521362 | -1.0343e-04 |
| 13 | 1.521362 | 1.521484 | 1.521423 | 2.5935e-04 |
| 14 | 1.521362 | 1.521423 | 1.521393 | 7.7956e-05 |
This method is a classic method that estimates a solution using the x-intercept of the line segment joining the endpoints of the function on the current interval. In this method, the root of the equation $f(x) = 0$ for the real variable $x$ is estimated when $f$ is a continuous function defined on a closed interval $[a, b]$ and where $f(a)$ and $f(b)$ have opposite signs. Moreover, $c_n$ is identified:
$$c_{n}={\frac {a_{n}f(b_{n})-b_{n}f(a_{n})}{f(b_{n})-f(a_{n})}}.$$def Regula_falsi(f, a, b, Nmax = 1000, TOL = 1e-4, Frame = True):
# endpoint values a, b,
# tolerance TOL,
# maximum iterations NMAX
if f(a)*f(b) >= 0:
print("Regula falsi method is inapplicable .")
return None
# let c_n be a point in (a_n, b_n)
Nmax=Nmax+1
an=np.zeros(Nmax, dtype=float)
bn=np.zeros(Nmax, dtype=float)
cn=np.zeros(Nmax, dtype=float)
fcn=np.zeros(Nmax, dtype=float)
# initial values
an[0]=a
bn[0]=b
for n in range(0,Nmax-1):
cn[n]= (an[n]*f(bn[n]) - bn[n]*f(an[n])) / (f(bn[n]) - f(an[n]))
fcn[n]=f(cn[n])
if f(an[n])*fcn[n] < 0:
an[n+1]=an[n]
bn[n+1]=cn[n]
elif f(bn[n])*fcn[n] < 0:
an[n+1]=cn[n]
bn[n+1]=bn[n]
else:
print("Bisection method fails.")
return None
if (abs(fcn[n]) < TOL):
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
n=n+1
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
Thus, for the above example, we have,
data = Regula_falsi(f, a, b, Nmax = 20, TOL = 1e-4)
#
Method_Plot(data, Name = 'Regula falsi', YL = [-1, .2], CL = 'RoyalBlue')
| an | bn | cn | fcn | |
|---|---|---|---|---|
| 0 | 1.000000 | 2.000000 | 1.333333 | -9.6296e-01 |
| 1 | 1.333333 | 2.000000 | 1.462687 | -3.3334e-01 |
| 2 | 1.462687 | 2.000000 | 1.504019 | -1.0182e-01 |
| 3 | 1.504019 | 2.000000 | 1.516331 | -2.9895e-02 |
| 4 | 1.516331 | 2.000000 | 1.519919 | -8.6751e-03 |
| 5 | 1.519919 | 2.000000 | 1.520957 | -2.5088e-03 |
| 6 | 1.520957 | 2.000000 | 1.521258 | -7.2482e-04 |
| 7 | 1.521258 | 2.000000 | 1.521344 | -2.0935e-04 |
| 8 | 1.521344 | 2.000000 | 1.521370 | -6.0461e-05 |
There have been also some improved versions of this method. For example, in the Illinois algorithm, $c_n$ is identifed as follows.
$$c_{n}=\dfrac {{\frac {1}{2}}f(b_{n})a_{n}-f(a_{n})b_{n}}{{\frac {1}{2}}f(b_{n})-f(a_{n})}$$def Regula_falsi_Illinois_alg(f, a, b, Nmax = 1000, TOL = 1e-4, Frame = True):
# endpoint values a, b,
# tolerance TOL,
# maximum iterations NMAX
if f(a)*f(b) >= 0:
print("Regula falsi (Illinois algorithm) method is inapplicable .")
return None
# let c_n be a point in (a_n, b_n)
Nmax=Nmax+1
an=np.zeros(Nmax, dtype=float)
bn=np.zeros(Nmax, dtype=float)
cn=np.zeros(Nmax, dtype=float)
fcn=np.zeros(Nmax, dtype=float)
# initial values
an[0]=a
bn[0]=b
for n in range(0,Nmax-1):
cn[n]= (0.5*an[n]*f(bn[n]) - bn[n]*f(an[n])) / (0.5*f(bn[n]) - f(an[n]))
fcn[n]=f(cn[n])
if f(an[n])*fcn[n] < 0:
an[n+1]=an[n]
bn[n+1]=cn[n]
elif f(bn[n])*fcn[n] < 0:
an[n+1]=cn[n]
bn[n+1]=bn[n]
else:
print("Bisection method fails.")
return None
if (abs(fcn[n]) < TOL):
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
n=n+1
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
data = Regula_falsi_Illinois_alg(f, a, b, Nmax = 20, TOL = 1e-4)
#
Method_Plot(data, Name = 'Regula falsi (Illinois algorithm)', YL = [-.15, .1], CL = 'Purple')
| an | bn | cn | fcn | |
|---|---|---|---|---|
| 0 | 1.000000 | 2.000000 | 1.500000 | -1.2500e-01 |
| 1 | 1.500000 | 2.000000 | 1.529412 | 4.8036e-02 |
| 2 | 1.500000 | 1.529412 | 1.524671 | 1.9614e-02 |
| 3 | 1.500000 | 1.524671 | 1.522877 | 8.9069e-03 |
| 4 | 1.500000 | 1.522877 | 1.522090 | 4.2212e-03 |
| 5 | 1.500000 | 1.522090 | 1.521723 | 2.0394e-03 |
| 6 | 1.500000 | 1.521723 | 1.521547 | 9.9423e-04 |
| 7 | 1.500000 | 1.521547 | 1.521462 | 4.8682e-04 |
| 8 | 1.500000 | 1.521462 | 1.521420 | 2.3888e-04 |
| 9 | 1.500000 | 1.521420 | 1.521399 | 1.1734e-04 |
| 10 | 1.500000 | 1.521399 | 1.521389 | 5.7666e-05 |
Interpolate Truncate and Project (ITP) method is an improved bisection method that can achive a superlinear convergence of the secant method. In worst-case, it matches the performance of the bisection method.
def ITP_Method(f, a, b, TOL = 1e-4, N0 = 1, K1 = 0.1, K2 = 2, Frame = True):
# endpoint values a, b,
# tolerance TOL,
# maximum iterations NMAX
Nmax = int(np.ceil(np.log2((b-a)/(2*TOL))) + N0)
# let c_n be a point in (a_n, b_n)
an=np.zeros(Nmax, dtype=float)
bn=np.zeros(Nmax, dtype=float)
cn=np.zeros(Nmax, dtype=float)
fcn=np.zeros(Nmax, dtype=float)
# initial values
an[0] = a
bn[0] = b
n = 0
while (bn[n]- an[n] > (2*TOL)):
mid = (an[n] + bn[n])/2
r = TOL * (2 ** (Nmax - n)) - (bn[n]- an[n])/2
xf = (an[n]*f(bn[n]) - bn[n]*f(an[n])) / (f(bn[n]) - f(an[n]))
sigma = np.sign(mid - xf)
delta = K1*((b-a)**K2)
#Truncation:
if delta <= np.abs(mid - xf):
xt = xf + delta * sigma
else:
xt = mid
# Projection
if np.abs(xt - mid) <= r:
cn[n] = xt
else:
cn[n] = mid - sigma * r
# Updating Interval:
fcn[n] = f(cn[n])
if fcn[n] > 0:
an[n+1]=an[n]
bn[n+1] = cn[n]
elif fcn[n] < 0:
an[n+1] = cn[n]
bn[n+1]=bn[n]
else:
an[n+1] = cn[n]
bn[n+1] = cn[n]
n += 1
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
Furthermore, for the above example, we have,
data = ITP_Method(f, a, b, TOL = 5e-4)
#
Method_Plot(data, Name = 'ITP Method', YL = [-.6, .6], CL = 'OrangeRed')
| an | bn | cn | fcn | |
|---|---|---|---|---|
| 0 | 1.000000 | 2.000000 | 1.433333 | -4.8863e-01 |
| 1 | 1.433333 | 2.000000 | 1.595020 | 4.6285e-01 |
| 2 | 1.433333 | 1.595020 | 1.514177 | -4.2576e-02 |
| 3 | 1.514177 | 1.595020 | 1.554599 | 2.0252e-01 |
| 4 | 1.514177 | 1.554599 | 1.534388 | 7.8091e-02 |
| 5 | 1.514177 | 1.534388 | 1.524282 | 1.7291e-02 |
| 6 | 1.514177 | 1.524282 | 1.519230 | -1.2759e-02 |
| 7 | 1.519230 | 1.524282 | 1.521756 | 2.2367e-03 |
| 8 | 1.519230 | 1.521756 | 1.520493 | -5.2684e-03 |
| 9 | 1.520493 | 1.521756 | 1.521124 | -1.5176e-03 |
| 10 | 1.521124 | 1.521756 | 0.000000 | 0.0000e+00 |